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We develop a relativistic lattice Boltzmann (LB) model, providing a more accurate description 
of dissipative phenomena in relativistic hydrodynamics than previously available with existing LB 
schemes. The procedure applies to the ultra-relativistic regime, in which the kinetic energy (tem- 
perature) far exceeds the rest mass energy, although the extension to massive particles and/or low 
temperatures is conceptually straightforward. In order to improve the description of dissipative ef- 
fects, the Maxwell- Jiittner distribution is expanded in a basis of orthonormal polynomials, so as to 
correctly recover the third order moment of the distribution function. In addition, a time dilatation 
is also applied, in order to preserve the compatibility of the scheme with a cartesian cubic lattice. 
To the purpose of comparing the present LB model with previous ones, the time transformation 
is also applied to a lattice model which recovers terms up to second order, namely up to energy- 
momentum tensor. The approach is validated through quantitative comparison between the second 
and third order schemes with BAMPS (the solution of the full relativistic Boltzmann equation), for 
moderately high viscosity and velocities, and also with previous LB models in the literature. Excel- 
lent agreement with BAMPS and more accurate results than previous relativistic lattice Boltzmann 
models are reported. 
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I. INTRODUCTION 



Relativistic hydrodynamics and kinetic theory play a 
major role in many forefronts of modern physics, from 
large-scale applications in astrophysics and cosmology, 
to microscale electron flows in graphene [H-Q, all the 
way down to quark-gluon plasmas 0-01 • to their 

strong non-linearity and, for the case of kinetic theory, 
high dimensionality as well, the corresponding equations 
are extremely challenging even for the most powerful nu- 
merical methods, let alone analytics. Recently, a promis- 
ing approach, based on a minimal form of relativistic 
Boltzmann equation, whose dynamics takes place in a 
fully discrete phase-space and time lattice, known as rel- 
ativistic lattice Boltzmann (RLB), has been proposed by 
Mendoza et al. 0-0- To date, the RLB has been ap- 
plied to the simulation of weakly and moderately rela- 
tivistic fluid dynamics, with Lorentz factors of 7 ^ 1.4, 
where 7=1/ \/l — j (? , c being the speed of light and 
V the speed of the fluid. This model reproduces correctly 
shock waves in quark-gluon plasmas, showing excellent 



agreement with the solution of the full Boltzmann equa- 
tion as obtained by Bouras et al. using BAM PS ( Boltz- 
mann Approach Multi-Parton Scattering) [13, [Hi- The 
RLB makes use of two distribution functions, the first 
one modeling the conservation of number of particles, 
and the second one, the momentum-energy conservation 
equation. The model was constructed by matching the 
first and second order moments of the discrete-velocity 
distribution function to those of the continuum equilib- 
rium distribution of a relativistic gas 
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In a subsequent work, Hupp et al.[l2| improved the 
scheme by extending the equilibrium distribution func- 
tion for the number of particles, in such a way as to 
include second order terms in the velocity of the fluid, 
thereby taming numerical instabilities for higher pressure 
gradients and velocities. However, the model was not 
able to reproduce the right velocity and pressure profiles 
for the Riemann problem in quark-gluon plasmas, for the 
case of large values of the ratio between the shear viscos- 
ity and entropy density, rj/s ^ 0.5, at moderate fluid 
speeds {v/c ^ 0.6). 

In order to set up a theoretical background for the lat- 
tice version of the relativistic Boltzmann equation, Ro- 
matschke et al. [31 developed a scheme for an ultrarel- 
ativistic gas based on the expansion on orthogonal poly- 
nomials of the Maxwell- Jiittner distribution [lj| and, by 
following a Gauss-type quadrature procedure, the dis- 
crete version of the distribution and the weighting func- 
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tions was calculated. This procedure was similar to the 
one used for the non-relativistic lattice Boltzmann model 
[isl . This relativistic model showed very good agree- 
ment with theoretical data, although it was not compat- 
ible with a lattice, thereby requiring linear interpolation 
in the free-streaming step. This implies the loss of some 
key properties of the standard lattice Boltzmann method, 
such as negative numerical diffusion and exact streaming. 

Very recently, Li et al. noticed that the equa- 

tion of conservation for the number of particles, recov- 
ered by the RLE model 0, Q , exhibits incorrect diffusive 
effects. Therefore, they proposed an improved version of 
RLE, using a multi-relaxation time collision operator in 
the Boltzmann equation, showing that this fixes the is- 
sue with the equation for the conservation of the number 
of particles. The generalized collision operator allows to 
tune independently the bulk and shear viscosities, yield- 
ing results for the Riemann problem closer to EAMPS 
[10| when the bulk viscosity is decreased. However, the 
third order moment of the equilibrium distribution still 
does not match its continuum counterpart and therefore 
the model still has problems to reproduce high rj/s ^ 0.5, 
for moderately high velocities, (3 = v/c = 0.6. Thus, 
while surely providing aji improvement on the original 
RLE model, the work |17| did not succeed in reproduc- 
ing the vanishing bulk viscosity, which is pertinent to the 
ultra-relativistic gas, while allowing the bulk viscosity to 
vary independently on the shear viscosity. 

Note that in the much more studied case of the lat- 
tice Boltzmann models for the non-relativistic fluids, the 
question of the choice of the lattice with higher-order 
symmetry requirements has only recently been solved, 
in the framework of the entropy-compliant constriction 
(Tsl [l9j . However, the lattices (space- filling discrete- 
velocity sets) found in that case are tailored to reproduce 
the moments of the non-relativistic Maxwell-Boltzmann 
distribution, and do not seem to be directly transferable 
to the present case of the relativistic (Maxwell- Jiittner) 
equilibrium distribution, which has fairly different sym- 
metries as compared to the non-relativistic Maxwell- 
Boltzmann distribution. Therefore, the extension of the 
previous LB models has to be considered anew. 

In this paper, we develop a new lattice Boltzmann 
model capable of reproducing the third order moment of 
the continuum equilibrium distribution, and still realiz- 
able on a cubic lattice. The model is based on a single dis- 
tribution function and satisfies conservation of both num- 
ber of particles and momentum-energy equations. The 
model is based on the single relaxation time collision op- 
erator proposed by Anderson and Witting |20| which 
is more appropriate for the ultra-relativistic regime than 
the Marie model used in the previous works. Thus, the 
proposed model offers significant improvement on previ- 
ous relativistic lattice Boltzmann models in two respects; 
(i) It captures the symmetry of the higher-order equilib- 
rium moments sufficiently to reproduce the dissipative 
relativistic hydrodynamics at the level of the Grad ap- 
proximation to the relativistic Boltzmann equation; (ii) 



It represents a genuine lattice Boltzmann discretization 
of space and time, with no need of any interpolation 
scheme, thereby avoiding the otherwise ubiquitous spu- 
rious dissipation. The new lattice Boltzmann model is 
shown to reproduce with very good accuracy the results 
of the shock-waves in quark-gluon plasmas, for moder- 
ately high velocities and high ratios ry/s. 

The paper is organized as follows: in Sec. [H] we de- 
scribe in detail the model and the way it is constructed; 
in Sec. IIIIl we implement simulations of the Riemann 
problem in order to validate our model and compare it 
with EAMPS and previous relativistic lattice Boltzmann 
models; finally, in Sec. IIVI we discuss the results and 
future work. 



II. MODEL DESCRIPTION 

A. Symmetries of the relativistic Boltzmann 
equation 

To build our model, we start from the relativistic 
Boltzmann equation for the probability distribution func- 
tion /: 



(1) 



where the local equilibrium is given by the Maxwell- 
Jiittner equilibrium distribution jl4{ . 



n = AcM-P^'U^/kBT) 



(2) 



In the above, ^ is a normalization constant, c the speed of 
light, and fc^ the Boltzmann constant. The 4-momentum 
vectors are denoted by = {E/c,p), and the macro- 
scopic 4-velocity by f/^ = {c,u)j{u), with u the three- 
dimensional velocity of the fiuid. Note that we have used 
the Anderson- Witting collision operator (rhs of Eq. 
([T|)). making our model compatible with the ultrarel- 
ativistic regime. Hereafter, we will use natural units, 
c = ks = 1, and work in the ultrarelativistic regime, 
i = mc^/kBT<^ 1. 

According to a standard procedure [H, [T^, [l^, we 
first expand the Maxwell- Jiittner distribution in the 
rest frame, /°'^ = Aexp(— p^/T), in an orthogonal 
basis. Since in the ultrarelativistic regime, /T ~ 
p/T, being p = \ 



fp' /T"^ + m? /T"^ ~ p/'J'j being p = \/ we can write 
the equilibrium distribution in spherical coordinates. 



gpe-P^^dpsm{9)ded(l) 



(3) 

where g is an arbitrary function of momentum. Follow- 
ing Romatschke |13j . we can expand the distribution in 
each coordinate separately, and subsequently, by using 
a Gauss quadrature, calculate the discrete values of the 
4-momentum vectors. Thus, the discrete equilibrium dis- 
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tribution can be written as, 

/^ = E«^^■fc(t^'^)^^(^O7^,(pz)Ffe(0o , (4) 

where the coefRcients aijk{U^^) are the projections of 
the distribution on the polynomials Pi{9i)TZj{pi)Fk{4'i), 
and the discrete 4-moinenta are denoted by = 
{pi,picos{4ii)sm{0i),pi sm{4ii)sm{0i),pi cos(6'/)). Conse- 
quently, the discrete form of the Boltzmann equation 
takes the form, 

TPi 

However, note that, in the streaming process on the right- 
hand-side of Eq.(l5]), the distribution moves at velocity 
pf/p^, which implies that the information travels (in a 
single time step) from each cell center to a position that 
belongs to the surface of a sphere of radius c St = 1. Fur- 
thermore, to represent correctly the third order moment 
of the equilibrium distribution, 

p'^^'' = j2frp?p?pt , (6) 
I 

the number of points needed on the surface of the unit 
sphere exceeds 6 and 12, which correspond to the first 
neighbors for a cubic and hexagonal closed packed (HCP) 
lattices, respectively. This implies that, in general, the 
4- vectors p^ /p^ cannot be embedded into a regular lat- 
tice, and therefore, an interpolation algorithm has to be 
used. By doing this, we are losing one of the most impor- 
tant features of lattice Boltzmann models, which is the 
exact streaming. Thus, within this spherical coordinate 
representation, the streaming process cannot take place 
on a regular lattice. 

B. Moment projection of the equilibrium 

In this work, we shall use a different approach to the 
quadrature representation. We first expand the equilib- 
rium distribution at rest, w{p'^) — f^^ ~ ^cxp(— p'^) 
by using Cartesian coordinates, unlike the spherical co- 
ordinate system used in Ref. [l^l, and choose the 4- 
momentum vectors such that they belong to the lattice 
(from now on and without loss of generality, we will use 
the notation p^/T p"). This procedure also avoids 
extra terms in the product, Pi{9i)TZj{pi)Fk{4>i) for the 
spherical case, which are not necessary if we only need to 
recover correctly the first three moments of the equilib- 
rium distribution. This considerably simplifies the dis- 
crete equilibrium distribution. 

By performing a Gram-Schmidt procedure with the 
weig ht w{p°), we construct a set of orthonormal polyno- 
mials. The orthonormal polynomials in cartesian coordi- 
nates up to third order, herefrom denoted by J^, where 
the index k runs from to 29, are shown in Table ID 



Order 


Polynomial Jk 


k 


0th 


1 





1st 


p"-2 p^' p» 


1, 2, 3, 4 


2nd 


(p"-6)p"+6 (p"-4)p^ (p"-4)p» 

2%/3 ' 2\/2 ' 2^2 
(pt'_4)p- _p02^px2_^2p„2 p02_3px2 

2/2 ' _ 4/2^ ' 4/6 
p^p" P^ P^ 
2-y2 ' 2v/2 


5, 6, 7 
8, 9, 10 
11, 12, 13 


3rd 


1 Cn" fi^2 o ((p"-10)p"+20)p" 
12 VP p Z, 

-^(/-6)(p°^-3p-),SE^:i^ 

((p"-10)p"+20)p!' (p°-6)p"=p" p==p«p= 
4/5 ' 4/3 ' 4/3 
(pO -6) (pO^ -p-2 -2p!'2 ) (_p02 +px2 +2p„2 ) 


14, 15 
16, 17 

18, 19, 20 
21, 22 
23, 24 

25, 26, 27 
28, 29 


8\/3 ' 8y/3 
py{-3p°^+3p^^+4py^) ((p0_l0)p0+20)p- 

24/2 ' 4/5 
(pO_6)p^p^ p" {p°^ -5p^^) (pO-ajp^p^ 
4/3 ' 8/30 ' 4/3 
(p°-6)pyp' p'(-p"2+p"^+4p''=) 
4^/3 ' 24,\/2 



TABLE 1: Polynomials Jk that are orthonormal on the weight 
function ■w{p''^) in Cartesian coordinates {x,y,z). 



Note that in this Table, the 4-momentum has the no- 
tation p^ = {p'^ ,p^ ,py ,p^). Since these polynomials are 
orthonormal, there are no normalization factors, and the 
Maxwell- Jiittner distribution can be approximated, up 
to third order in the momentum space, by an expansion 
as simple as 

29 

r'i:^^«;(pO)afc(T,[/'^)Jfe(p^) , (7) 

fe=0 

where the projections are calculated by, 

a. = I n-Jkipn^ . (8) 

Since the Anderson- Witting model is only compatible 
with the Landau-Lifshitz decomposition [l^.l20l| , we must 
calculate the energy density of the fluid by solving the 
eigenvalue problem, 

T°'^Up eU"' , (9) 
e being the energy density of the fluid, and 

T-'= JfP^p'^ . (10) 

the momentum-energy tensor. For the particle density, 
we use the relation, 

n^U^Jfp'^^ , (11) 

and, by using the equation of state, e = 3nT, we can 
calculate the temperature of the fluid. 



C. Discrete-velocity representation of the 
quadratures 

Note that the above derivation using Cartesian coordi- 
nates still refers to the continuous 4-momenta. In order 



4 



to discrctize the above moment projection of the equihb- 
rium distribution, we must choose a set of 4-momentum 
vectors that satisfies the very same orthonormahty con- 
ditions, namely: 

^ i 

(12) 

while, at the same time, p'^/p'^ corresponds to lattice 
points. Here, we choose to work with a cubic lattice, al- 
though the procedure described here also applies to other 
ones, e.g. HCP lattice. 

Since, due to its nature, p^/p'^ leads to velocity vectors 
which belong to a sphere of radius c in the space com- 
ponents, using the procedure in Ref. [l3| will generally 
result in off-site lattice points. For this reason, we opt for 
another quadrature based on this orthonormahty condi- 
tion, and impose that the distribution function at rest 
frame should satisfy the moments of the equilibrium dis- 
tribution, up to 6-th order. This is made to ensure that 
the 5-th order moment of the equilibrium distribution is 
recovered (at least at very low fluid velocities), which, in 
the context of the Grad theory for the Anderson- Witting 
model [l3| , is a requirement for the correct calculation of 
the transport coefficients, namely the shear and bulk vis- 
cosities and thermal conductivity. The condition for the 
6-th order moment, is to choose from the multiple lattice 
solutions, the one that presents the highest symmetry to 
model the Maxwell- Jiittner distribution. In order to use 
general features of classical lattice Boltzmann models, 
like bounce-back boundary conditions to impose zero ve- 
locity on solid walls, wc will also require that the weights 
Wi corresponding to the discrete 4-momentum vectors pf 
have the same values as the ones corresponding to — pj^ 
(latin indices run over spatial components). 

In order to generate on-site lattice points, let us first 
analyse the relativistic Boltzmann equation, which can 
be written as, 

pOdJ+p-daf=-^^if-n) , (13) 
T 

and in the ultrarelativistic regime, 

pO(dtf+v'^daf) = -^{f-n) . (14) 
T 

where w° are the components of the microscopic velocity. 
These microscopic velocities have the same magnitude 
but, in general, different directions. Dividing both sides 
of Eq. lfTB)) by p", we obtain 

dtf + v-dj^~^i^{f-n) . (15) 

Tp^ 

In other words, in the ultra- relativistic regime, the rela- 
tivistic Boltzmann equation can be cast into a form where 
the time derivative and the propagation term become 
the same as in the non-relativistic case, at the price of 
an additional dependence on p^ in the relaxation term. 




FIG. 1: Directions of the velocity vectors "di to recover up to 
the third order moment of the Maxwell- Jiittner distribution. 
The radius of the sphere is i? = \/41. The points represent 
lattice sites belonging to the sphere surface. 

However, since this newly acquired dependence remains 
local, we shall be able to find a discrete- velocity quadra- 
ture which also allows for a lattice Boltzmann-type dis- 
cretization in time and space without any interpolation. 
Indeed, in a cubic cell of length Sx — I there arc only 6 
neighbors, which are not sufficient to satisfy the orthog- 
onality conditions and the third order moment of the 
equilibrium distribution. However, by multiplying this 
equation by a constant R at both sides, and perform- 
ing a time transformation (dilatation), St — )■ R6t' and 
T Rt' , we obtain 

dt'f + rdaf = -^{f-n) , (16) 

where we have defined "(9° = Rv'^. Due to this transforma- 
tion, the 4-momentum vectors are reconstructed through 
the relation 

p'^^p°{lJ/R) , (17) 

At this stage, we can choose the radius of the sphere such 
that the lattice points that belong to the surface of the 
sphere and the cubic lattice exhibit enough symmetries 
to satisfy both conditions. This is equivalent to solving 
the Diophantine equation, 

nl + nl + nl^R' , (18) 

where n^, Uy, and are integer numbers, being t? = 
(jixjUyjUz). Thus, wc can determine the components of 
the discrete version of the velocities which are needed 
for the streaming term in the Boltzmann equation. Ihs of 
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Eq. p6)) . However, on the rhs of this equation, and for 
the calculation of the discrete 4-monientuni vectors via 
Eq. (fT7|) . we also need to know the discrete values of p". 
The 4- vector is needed to compute the orthonormality 
conditions given by Eq. (|12p and the moments of the 
equilibrium distribution. 

Due to the fact that p° is the magnitude of the 4- 
momentum, p° = ^Jp^^p^ , in 3 + 1-dimensional space- 
time, it is natural to assume that its discrete values can 
be calculated by using the weight function in spherical 
coordinates, w{p) = 47rAp^ exp(— p), where the angular 
components have been integrated out, and using the zeros 
of its respective orthonormal polynomial of fourth order 
(this is because we are interested in an expansion up to 
third order, so we need one more order to calculate the 
zeros). This fourth order polynomial is given by: 



w, > 



(22f) 



1 



24 



[120 + p(-240 + p[120 +{p- 20)p])] 



(19) 

To summarize, in order to calculate the discrete p^ and 
their respective Wi , we first fix R and solve the equations 



7^W(p) = 



(20a) 



(20b) 



to obtain the solutions for rix, riy, n^, and p. With these 
values, we build the discrete 4- vectors 



Pfm ^ PU^^'^x,m/R,ny,jn/R,n^^rn/Ii) 



(21) 



where Z = 1, 4 denotes the four zeros of the polynomial 
7?.*-^^(p), and m = 0,...,A^ the triplets {nx,ny,nz)m that 
satisfy the Diophantine equation, assuming that M is the 
number of solutions. Here, for simplicity, we regroup the 
pair of indexes im to i, so that we can label the discrete 
4-momentums as pf , where i = 1, A/" with Af = Ax Ai. 
Next, we replace these values into the equations. 



J w{p\Ji{pn.h{p^)^ - y^MDMpf) = Sik , 

(22a) 



w{p )p'"p p p = 2_^w.ip'^p.^p^p. 



J «;(pO)p>Vp^p''-i: = Y^w,p1p^p^p: 



'p1 



(22b) 



(22c) 



and look for any solution for Wi that fulfills the above 
relations. Should none be found, we repeat the procedure 
with a different value of R. By performing this iteration 
process, we found that R = ^/41 is sufficient to recover 
up to the third order moment of the Maxwell- Jiittner 
distribution, and up to sixth order of this distribution in 
the Lorentz rest frame. 

The corresponding discrete velocity vectors 79„j are: 
(±6,±2,±1), (±6,±1,±2), (±2,±6,±1), (±1,±6,±2), 
(±1,±2,±6), (±2,±1,±6), (±5,0,±4), (±5,±4,0), 
(0,±5,±4), (±4, ±5,0), (0,±4,±5), (±4,0, ±5), 
(±4, ±3, ±4), (±3, ±4, ±4), and (±4, ±4, ±3); with 
the values for p° ~ 0.743, 2.572, 5.731, and 10.95. 
Consequently, this gives a total of 4-momentum vectors 
TV = 384. However, the last condition in Eq. ((^ 
allows some weights to become zero. Therefore, in our 
iteration procedure, we have taken the minimal number 
of 4-momentum vectors p^, by imposing the maximum 
number of Wi to be zero. For this reason, there are 
only 128 vectors pf needed to fulfill the conditions in 
Eq. (|22p . In principle, all the velocity vectors t?„i are 
needed, but only some of the combinations with pj* are 

required. The detailed list of the p^, and p^, and 
their respective discrete weight functions Wi are given in 
the Supplementary Material [2]| . 

In Fig. [1] we report the configuration of the veloc- 
ity vectors to achieve the third order moment of the 
Maxwell- Jiittner distribution function. The points cor- 
respond to lattice nodes of a cubic lattice that, at the 
same time, belong to the surface of the respective sphere 
of radius R = \/41. The relatively large number of dis- 
crete velocities should not come as a surprise; in the case 
of non-relativistic lattice Boltzmann, the number of dis- 
crete velocities also becomes high (at least 41 for achiev- 
ing complete Galilean invariance in the non-thermal case 
and 125 in the thermal case, see [H, [I!]). Note that the 
specified values of p'^ play the same role in defining the 
quadrature as the reference temperature (energy) in the 
non-relativistic case [3 [3 ■ 

Finally, we can write the discrete version of the equi- 
librium distribution up to third order. 



29 

E 

n=0 



(23) 



which is shown in details in Appendix iBl Eq. (jB2[) . Note 
that this distribution function recovers the first three mo- 
ments of the Maxwell- Jiittner distribution in the ultra- 
relativistic regime. 



// Q \ n V (J X ^ R '^^P u y 

w{p )p'"p p p P^P ^ 2_^w^p'^p, 



cr A 7 /3 
Pt Pt PiPi > 



Wi 



(if P- = -Pj ) , 



(22d) 
(22e) 



,Q 128 



cq u 

Pi 



ry^p' 



,3 128 

^ i=l 



cq II V 
PiPi 



(24) 



(25) 
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rwp'^^Y.f^''p'p'^p' = p"^' , (26) 



where 



(27) 



(28) 



being the number of particles 4-fiow and the energy- 
momentum tensor, respectively, and 

(29) 

with n ~ 2T^ . However, the extension to the case of 
massive particles is straightforward, by changing the co- 
efficients, a„, in Eqs. ((8)) and (|23)) . 



D. Discrete relativistic Boltzmann equation 

In the model of Anderson- Witting for the collision 
operator, the relativistic Boltzmann equation takes the 
form given by Eq. ([T}, 



(30) 



This collision operator is compatible with the Landau- 
Lifshitz decomposition [l^ . which implies fulfillment of 
the following relations 

(31a) 

U^T^-" = U^J fp^p"^ = U^T^- = J np^p"^ , 

^ . . 

Here, the subscript E denotes the quantities calculated 

with the equilibrium distribution. Therefore, upon inte- 
grating Eq. pop in momentum space, we obtain 



(32) 



which is the conservation of the number of particles 4- 
fiow. By multiplying by p"^ and integrating, we obtain 
the conservation of the momentum energy tensor 



d^Tt"" ^ 



(33) 



In order to calculate the transport coefficients, we need 
the third order moment, so that, upon multiplying 
Eq. (1201) by p'^p'^, we obtain 



(34) 



and by using a Maxwellian iteration method [l^ , 



E 



(35) 



Note that we need at least the third order moment of 
the equilibrium distribution, P^'^^ , to compute the dis- 
sipation coefficients (namely, bulk and shear viscosities 
and heat conductivity). This requirement is fulfilled in 
our discrete and continuum expansions of the equilib- 
rium distribution via Eqs. p4l). (l25l) . (pHl) . However, to 
recover full dissipation, we would also need to recover the 
third moment of the non-equilibrium distribution, which 
according to the 14 moments Grad's theory , can be writ- 
ten as. 



(36) 



where ba and da\ are coefficients that carry the infor- 
mation on the transport coefficients [l3|- Note that we 
need to recover terms up to the fifth order of the equi- 
librium distribution. In principle, this could be done by 
the procedure described on this paper, but the resulting 
value for R could be unpractically large. Nevertheless, 
at low velocities, C/^ ^ (1,0,0,0), the Maxwell- Jiittner 
distribution can be approximated by the weight func- 
tion w{p'^), and in analogy to the discrete case, by Wi, 
and the fourth and fifth order are recovered via Eq. (|22p . 
As a result, at relatively low velocities, we expect the 
non-equilibrium third order tensor to be also fulfilled. 
Therefore, the transport coefficients for an ultrarelativis- 
tic gas, i.e. = for the bulk viscosity, 77 = (2/3)Pt for 
the shear viscosity, and A = (4/5T)Pr for the thermal 
conductivity, also apply to our model. 

To discretize the relativistic Boltzmann equation, we 
first implement the time transformation described in the 
previous section and integrate in time Eq. pOI) between 
t' and t' + St'. This yields: 

t'p^ 

(37) 

By changing p^^ ^ p^^ , f ^ fi and -d"- Sf, we obtain 

t'p" 

(38) 

This relativistic lattice Boltzmann equation presents an 
exact streaming at the left hand side, and the collision 
operator at the right hand side looks exactly like its con- 
tinuum version. Therefore, the conservation laws for the 
number of particles density 4-flow, and the momentum- 
energy tensor, are also fulfilled, as long as they are ob- 
tained by using the Landau-Lifshitz decomposition. This 
means that, first, we need to calculate the momentum- 
energy tensor. 



rpaf: 



128 

E 

i=l 



Mp^ 



(39) 
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and with this tensor, we solve the eigenvalue problem, 

T^'PUp = T'^'^Uf} = eiJ" , (40) 

obtaining the energy density e and the 4- vectors 
Subsequently, the particle density can be calculated by 



128 



(41) 



The temperature T is obtained by using the equation of 
state for the ultrarelativistic gas, e ~ 3nT. The trans- 
port coefficients are the same as in the continuum case, 
with the lattice correction resulting from second order 
Taylor expansion of the streaming term. All factored in, 
the coefficients take the following expression = 0, r] = 
(2/3)P(t' - St'/2), and A = (4/5r)F(T' - St'/2). Note 
that reverting back the time transformation, we can write 
the transport coefficients as = {2/3)P{T — 6t/2)/R, and 
A ^ (4/5T)P(t - 5t/2)/R. 

Summarizing, the present model does not present spu- 
rious dissipation in the number of particle conservation 
equation, in contrast to previous RLB schemes 0, H, [l^l ; 
and also improves the dissipative terms given by the 
multi-relaxation time scheme [TtI . In addition, it realizes 
the expansion of the Maxwell- Jiittner distribution on a 
cubic lattice, in contrast to Ref. [l^. We can also con- 
struct a relativistic lattice Boltzmann model that recov- 
ers only up to second order (momentum-energy tensor), 
to compare with the third order model and determine 
the influence of the third order moment in the expan- 
sion. Details of the second order model can be found in 
Appendix \^ 
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FIG. 2: Velocity (top) and pressure (bottom) profiles as func- 
tion of the ^-coordinate for the case of a Shockwave in quark- 
gluon plasma, with rj/s = 0.1. 



III. NUMERICAL VALIDATION 

In order to validate our model, we solve the Riemann 
problem for a quark-gluon plasma and compare the re- 
sults with BAMPS and two previous relativistic Boltz- 
mann models. The first one, proposed by Mendoza et al. 
0, iS] and later improved by Hupp et al. [T^], which we 
will denote simply by RLB, and the second one, which is 
a recent extension of the RLB developed by Li et al. 
to include multi-relaxation time, which we will denote by 
MRT RLB. BAMPS was developed by Xu and Greiner 
and applied to the Riemann problem in quark-gluon 
plasma by Bouras et al. [ll|. Since BAMPS solves the 
full relativistic Boltzmann equation, we take its result as 
a reference to access the accuracy of our model. However, 
we keep in mind that BAMPS also produces approxi- 
mate solutions. The present model is hereafter denoted 
by RLBD (RLB with Dissipation). 

For small ratios rj/s, where s is the entropy density, 
RLB and MRT RLB reproduced BAMPS resuhs to a 
satisfactory degree of accuracy. However, for higher 
■q/s > 0.1 and moderately fast fluids, 7 ~ 1.3, RLB 
failed to reproduce the velocity and pressure proflles 



[1^ . MRT RLB yielded good agreement with the re- 
sults at rj/s = 0.1, but presented notable discrepancies 
for fj/s = 0.5. The failure of both RLB and MRT RLB 
to solve the Riemann problem for high viscous fluids can 
be ascribed to their inabilityto recover the third order 
moment of the distribution |T2| . [TtI . 

In this section, we will study the case of high rj/s > 0.1 
in a regime of moderate velocities. We perform the sim- 
ulations on a lattice with 1 x 1 x 1600 cells, only half of 
which are represented in our domain owing to symme- 
try condition (the other half is a mirror, in order to use 
periodic boundary conditions for simplicity). Therefore, 
our simulation consists of 1 x 1 x 800 lattice sites, with 
5x ^ 0.008/m and 5t = ^41 0.008/m/c for RLBD third 
order, and 5t = 0.024 fm/c for RLBD second order. 

The initial conditions for the pressure are Pq = 
5.43^ and Pi = 0.339^. In numerical units, they 
correspond to 1.0 and 0.062, respectively. The initial 
temperature z > is Ti = 200MeF (in numerical units 
0.5), and Tq = AOOMeV for z < 0, which corresponds 
to 1.0 in numerical units. The entropy density s is cal- 
culated according to the relation, s = An — nhi{n/n'^^), 
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FIG. 3: Velocity (top) and pressure (bottom) profiles as func- 
tion of the a-coordinate for the case of a Shockwave in quark- 
gluon plasma, with t)! s = 0.5. 



where is the density calculated with the equilibrium 
distribution, = (IqT^ /ii'^, with do = 16 being the 
degeneracy of the gluons. 

The velocity and pressure profiles at t = 3.2^ with 
viscosity-entropy density ratios of rj/ s = 0.1, are shown 
in Fig. [5] In this figure, we compare the results with 
BAMPS and RLB, where we can sec that RLB presents a 
discontinuity at z = 0, while both second order and third 
order RLBD get closer to the BAMPS solution. Since the 
only difference between second and third order RLBD is 
the third order moment of the distribution, we conclude 
that at relatively low r//s, the third order does not play a 
crucial role neither in the conservative dynamics nor dis- 
sipative dynamics of the system. However, note that at 
z 3fm, the third order model provides an outstanding 
fit of the numerical results by BAMPS. 

On the other hand, by increasing the ratio ?//s, we see 
from Fig. |3] that, while RLB gets worse and the second 
order RLBD fixes the discrepancy only in part, the 3rd 
order RLBD improves significantly the accuracy of the 
velocity and pressure profiles. 

In Fig. 131 we also compare the results obtained with 
MRT RLB and BAMPS, for r]/s = 0.5. Here, we see 




FIG. 4: Velocity profile as function of the z-coordinate for the 
case of a Shockwave in quark-gluon plasma, with rj/s = 0.5, by 
increasing the reference numerical temperature in the lattice, 
leading to a smaller relaxation time r. 



that there is again an improvement, including the at- 
tainment of the right value of the maximum velocity (at 
z ^ 1.5/m). In the pressure profile, RLBD gets closer to 
BAMPS than MRT RLB in the region of the discontinu- 
ity in the initial condition (z ^ 0). 

Note that there is a staircase shape in the results of 
RLBD for rj/s = 0.5 in Figs. [3l This is due to the large 
values taken by the single relaxation time in order to 
achieve such shear viscosity-entropy density ratios, r 
20 — 40 (in numerical units) , which is beyond the hydro- 
dynamic approximation and therefore higher order mo- 
ments (fourth and higher orders) of the distribution func- 
tion would be required, which is not fulfilled in our RLBD 
model. In order to prove this statement, we have per- 
formed separate simulations, see Fig.|4l where we observe 
that by increasing the value of the reference temperature 
of the lattice (typically set at T = 1), so as to achieve 
the same shear viscosity, 77 = (2/3)rtT(r — 1/2) /i?, the 
value of T decreases and the staircase disappears. In par- 
ticular, for T > 2.5, the results get closer to the ones 
with BAMPS, and become independent of the reference 
temperature. Unfortunately, due the discretization pro- 
cedure used to develop this model, whenever the refer- 
ence temperature T > 4 the model becomes unstable, 
mostly likely because the expanded equilibrium distribu- 
tion function takes negative values. 



IV. CONCLUSIONS AND DISCUSSIONS 

We have introduced a new relativistic lattice Boltz- 
mann model with improved dissipation, as compared to 
RLB and MRT RLB. To this purpose, we have per- 
formed an expansion of the Maxwell- Jiittner distribu- 
tion onto an orthonormal basis of polynomials in the 
4-momentum space. In addition, in order to make the 



9 



model compatible with a regular cubic lattice, we have 
performed the expansion in cartesian coordinates and ap- 
plied a time transformation, such that particles travel 
just the distance necessary to reach lattice nodes, always 
at the speed of light. The time transformation generates 
a sphere of radius R which intersects the cubic lattice, 
the intersection points being lattice nodes by construc- 
tion. In addition, wc have reproduced up to second order 
moment of the equilibrium distribution, and up to third 
order moment, finding i? = 3 and R = -\/4T for second 
and third order moment compatibility, respectively. 

The discrete energy component of the 4-momcntum, 
has been calculated by using Gaussian quadrature, 
the nodes corresponding to the zeros of the next order 
polynomial. With this configuration, we need 90 vectors 
for recovering second order and 384 for the third order 
moment case. However, only 66 and 128, respectively, 
are actually needed to calculate the moments correctly. 

In order to validate the model, we have compared our 
results with BAMPS, as well as previous RLB models. 
We have found that for ry/s = 0.1, our model accurately 
describes the Riemann problem in quark-gluon plasma, 
including the expansion up to second order. However, for 
the case of rj/s = 0.5, the second order model, although 
better than RLB, is less accurate than both MRT RLB 
and the third order model. The third order model yields 
better results than the previous RLB, but it develops a 
staircase shape as a consequence of the large value of the 
single relaxation time, which lies beyond the hydrody- 
namic regime. We have shown that the staircase pathol- 
ogy can be tamed by increasing the reference tempera- 
ture in the model. Nevertheless, increasing the reference 
temperature beyond T = 4 hits against stability limits of 
the model. 

Wc may envisage that a multi-relaxation time exten- 
sion of the present model would further improve the ac- 
curacy of the results. A similar improvement may be 
anticipated by implementing higher order expansions of 
the equilibrium distribution. However, since the trans- 
port coefficients depend on the collision operator, their 
calculation within a multi-relaxation time model becomes 
increasingly involved. On the other hand, by performing 
expansions to include higher order moments, the value of 
R might become unpractically large, with several ensu- 
ing discretization issues. Notwithstanding such potential 
difficulties, these extensions are surely worth being ana- 
lyzed in depth for the future. 
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FIG. 5: Directions of the velocity vectors §i to recover up 
to the second order moment of the Maxwell- Jiittner distribu- 
tion, namely the momentum-energy tensor. The radius of the 
sphere is i? = 3. The points represent lattice sites belonging 
to the surface of the sphere. 



Appendix A: Second order relativistic lattice 
Boltzmann model 



To construct the second order lattice Boltzmann 
model, we use the procedure described in this paper. We 
have obtained that i? = 3 presents enough symmetries 
to fulfill the conditions in Eqs. ([2^ . and the velocity 
vectors are given by, (±3,0,0), (0,±3,0), (0,0, ±3), 
(±2,±1,±2), (±1,±2,±2), and (±2,±2,±1). The val- 
ues for the discrete p° come from the solution of the equa- 
tion. 



1 



7e(3) = _L/(/ _ 6)2 



12 







(Al) 



instead of TZ^^^ for the case of the third order expan- 
sion. This gives the values p° ~ 0.936, 3.305, and 
7.759. The discrete 4-momentum vectors are con- 
structed with Eqs. (|17p . and (|2ip . and they are in total, 
A/" = 3 X 30 = 90. However, as in the third order ex- 
pansion, we have retained the minimal amount, out of 
90, that are necessary to recover the second order mo- 
ment, by imposing the maximum number of Wi to be 
zero. This gives only 66 4-momentum vectors. The value 
of the weight functions for every momentum vector and 
the relation with the 30 directions are given in the Sup- 
plementary Material [2l|. In Fig. [5] we report the spatial 
configuration of the vectors 'di. 

The discrete version of the relativistic Boltzmann equa- 
tion, Eq. (|38p . still applies and the discrete equilibrium 
distribution function is written in detail in Appendix |B1 
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Eq. (jBip . However, due to the fact that the third order 
moment is not satisfied, an analytical theory to calculate 
the transport coefficients would be very complicated and 
goes beyond the scope of this work. Therefore, we have 
calculated numerically only the shear viscosity, by match- 
ing the results for low velocity with the third order mo- 
ment model. This, in order to compare the results of both 
expansions with other models in the literature. This gives 
a shear viscosity 772nd ^ {^n)P{''' ~ St/2)/R. We could, 
in principle, calculate the third order moment associated 
with the equilibrium distribution given by Eq. (|B1[) , and, 
by applying the Grad method, compute the other trans- 
port coefficients. However, this procedure would need 
to be performed entirely numerically, since the weights 
Wi and 4-momentum vectors are only known numeri- 
cally. Since the main purpose of this paper is to improve 



the description of dissipative effects by performing the 
third order expansion and place it on a cubic lattice, we 
arc not interested in the bulk viscosity and the thermal 
conductivity for this case, and leave this task for future 
work. 



Appendix B: Equilibrium Distribution Functions 

The equilibrium distribution function capable to re- 
cover the first and second order moments of the equilib- 
rium distribution is calculated by using up to the second 
order polynomials in Eq. ([7]), namely the 14 polynomials 
Jk with fc = 0, 13, obtaining 



4T 



2TU" + 1) + 2p'^{T{T{U°{p^U'' +p^^Uy +ptU'' - 4U") + 1) 



+ pf ( -c/"^ + + 2uy'' + 1 ) + 2pfuy{p', u' - 4c/") + m"iu" - ptu') - 2 



+ 2r(5(pfc/" -t-pft/^ +pfc/") - St/") -t- 12 



(Bl) 



For the case of the third order moment expansion, we (fc = 0, ...,29). This leads to the following expressions: 
repeat the same procedure, using all the polynomials 
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nwi 



J i 1 rtrr 



12T 

2 



pt'{TU" - 1) ( ( - 3 ( C/"^^ + C/y^ + 1 ) ) - 2rC/" + 1 



- 3r(pf [7=^ + pf + pf - 14[/°) - 15^ - 3p^ (t^ (^U"^ (pf + pf - 24) - (^pf (2^"^ + U"^ + l) 

+ 2p^u^{pyuy +p'-u') +p\ (p'fu''^ + 2pyuy^+pf + 2p^t/^;7^~) - 12^ + uu^^p^w" +pyuy +ptu') 

- 2(pf + pfuy+piu')^ + (^f + 2c/-"' + + 1) + 2p^u^pyuy +p^u' - nu°) 

+ pf (-[/°^ + ?7^^ + 2[/«2 + 1) + '^P^Uy{ptU'- ~ 11U°) - 22plU°U' + 56U°^ - 14^ 
+ 2r(6(p;^[/^ + p^f/^ +pf [/'') - 25[/°) + 20^ +T[pfT'^U^ (-3[/"^ + 4J7^2 + 3?7'^2 + s) 
+ pf T (3 ([/"^ - 2C/-^2 - [/^^ - 1) i-pyruy + 6TU" - 7) +pITU'- (-U°^ + 4t/^2 ^ ^ -^^^ 
+ 3pf C/^ f T f T (pf ^ f-C/"^ + [/^^ + 2[/^^ + 1) + 2pyuy{plU'- - 6J7") - UplU^U"- + 24J7°^ - 4 



+ lApyuy + 14pf C/^ - 481/° j + 30 j + pfr'^uy (-3[/°^ + 3C/^^ + 4[/'^^ + 3) 
+ pf T (pfTt/^ (-[/"^ + J7^^ + 4J7f2 + 1) + 3(6rC/" - 7) (?7"^ - [/^^ - 2[/^2 - l)) 
+ Qpfuy (t (2T i^-iplU^U^ + 6C/"^ " 1) + ~ 24L/") + is) + 6p^t/^ (^T"^ {qU^ 

- 24U° (t^ (2U"^ + 5) ) - 30 (T^ - 2) 



1 ) - 24TU" + 15) 

m 
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